##############################################################
#Written by Collin Storlie
#This script is to bring in data from an Access Database
#and perform a Quantile Regression between two variables

#list the libraries needed

necessary=c("RODBC","quantreg")

#check if library is installed

installed = necessary %in% installed.packages()

#if library is not installed, install it

if (length(necessary[!installed]) >=1) install.packages(necessary[!installed], dep = T)

#load the libraries

for (lib in necessary) library(lib,character.only=T)

#define the input data directory

working = "C:/Users/Collin/Documents/PhD/R/Working/"

#define the output folder

outfolder = "C:/Users/Collin/Documents/PhD/R/R Outputs/"

#set the working directory where the asc files are stored

setwd(working)

#create a connection to Access Database

soildatabase <- odbcConnectAccess("C:/Users/Collin/Documents/PhD/Databases/MVD/MVD_Microclimate.mdb")

#Extract All Soil Moisture Data

t.query1 = "SELECT Data.soil_moisture, Data.date_time, georeference.point_ID FROM georeference INNER JOIN (Missions INNER JOIN Data ON Missions.mission_id = Data.mission_id) ON georeference.georef_ID = Missions.georef_id;"
t.soilmoisturedata = sqlQuery(soildatabase,t.query1)

#Extract All Daily Rainfall Data

t.query2 = "SELECT georeference.point_ID, Daily_Rainfall_Formatted.date, Daily_Rainfall_Formatted.daily_rainfall FROM Daily_Rainfall_Formatted INNER JOIN georeference ON Daily_Rainfall_Formatted.point_ID = georeference.point_ID;"
t.rainfalldata = sqlQuery(soildatabase,t.query2)

#Extract All Soil Moisture Linear Regression Stats

t.query3 ="SELECT [Soil Moisture Calibration Values].point_id, [Soil Moisture Calibration Values].p, [Soil Moisture Calibration Values].rsquare FROM [Soil Moisture Calibration Values];"
t.linregressstats = sqlQuery(soildatabase,t.query3)

#close the connection to database
odbcClose(soildatabase)

#Append Year and Month to All Data

year =(substr(t.soilmoisturedata$date_time,0,4))
month = (substr(t.soilmoisturedata$date_time,6,7))
t.soilmoisturedata2 = cbind(t.soilmoisturedata, year=year, month=month)
year2 =(substr(t.rainfalldata$date,0,4))
month2 = (substr(t.rainfalldata$date,6,7))
t.rainfalldata2 = cbind(t.rainfalldata, year=year2, month=month2)

#Produce Summary Table of Monthly Soil Moisture & Sum of Rainfall Data

daysrain = function(x) {
  numdays = which(x>2)
  lengthnumdays=length(numdays)
  return(lengthnumdays)
}

# daysrain = function(x) {return(length(which(x>2)))}

monthlysoil = aggregate(t.soilmoisturedata2$soilmoisturecalibrated,by=list(point_ID=t.soilmoisturedata2$point_ID, year=t.soilmoisturedata2$year, month=t.soilmoisturedata2$month),FUN=mean)
sumrain = aggregate(t.rainfalldata2$daily_rainfall, by=list(point_ID=t.rainfalldata2$point_ID, year=t.rainfalldata2$year, month=t.rainfalldata2$month),FUN=sum)
daysrainpermonth = aggregate(t.rainfalldata2$daily_rainfall, by=list(point_ID=t.rainfalldata2$point_ID, year=t.rainfalldata2$year, month=t.rainfalldata2$month), FUN=daysrain)


#Combine Summary Table of Monthly Soil Moisture & Sum of Rainfall

outdata = merge(monthlysoil, sumrain, by.x = c("point_ID","month","year"), by.y = c("point_ID","month","year"))

#Loop to create plots for each site with regression line through lower quartile and a summary table of regression coefficients for each site

site.list = unique(outdata$point_ID)

models.coefs = data.frame(point_ID=NA, slope=NA, intercept=NA)

for(site in site.list) {

tdata = subset(outdata,point_ID==site)

quantregress = rq(tdata$x.x~tdata$x.y,tau=0.1)

models.coefs = rbind(models.coefs,data.frame(point_ID=site, intercept=coefficients(quantregress)[1], slope=coefficients(quantregress)[2]))

png(paste(outfolder,site,".png"))

plot(tdata$x.y,tdata$x.x,main=site,xlab="Rainfall", ylab="Soil Moisture")

abline(quantregress)

dev.off()

}

#Perform Quantile Regression at 5%

quantregress = rq(outdata$x.x~outdata$x.y,tau=0.05)

plot(outdata$x.x~outdata$x.y)
abline(quantregress)